# Load libraries
library(data.table)
library(sctransform)
library(Seurat)
library(tidyverse)
library(miQC)
library(SeuratWrappers)
library(fishpond)
library(Matrix)
library(patchwork)
library(alevinQC)
library(harmony)
library(scDblFinder)
# Set global options for Seurat v5 objects
options(Seurat.object.assay.version = 'v5')2 IFNg: Normalization, doublet discrimination, integration, clustering
2.1 Set up Seurat workspace
2.2 Load previously saved merged object
merged.18279 <- readRDS("IFNg_scRNA_Part1.rds")2.3 Normalize and scale data
Regress out mitochondrial contribution
merged.18279 <- NormalizeData(merged.18279)
merged.18279 <- FindVariableFeatures(merged.18279,
assay = "RNA",
layer = "data",
selection.method = "vst",
nfeatures = 5000)
merged.18279 <- ScaleData(merged.18279, vars.to.regress = "percent.mt")2.4 Run PCA
merged.18279 <- RunPCA(merged.18279, npcs = 200)PC_ 1
Positive: FABP5, PIM3, TNFRSF9, CD82, TNFRSF4, GNG4, PGAM1, PRDX1, ATP1B3, RAN
DUSP4, TNFRSF18, RPS17, NME1-NME2, C1QBP, MIR155HG, NDFIP1, NPM1, RPL7, RANBP1
ID2, SRM, HSP90AB1, HSPE1, PARK7, BACH2, PDCD1, HSPD1, REL, MAP2K3
Negative: IL32, CCL5, PLAAT4, LTB, TMSB4X, GZMA, CD74, PTPRC, PTPRCAP, CORO1A
IFITM1, MYL12A, CYBA, CYTIP, ANXA1, IFI6, GIMAP7, CD8A, RNF213, EVL
TMSB10, CD8B, LIMD2, SAMHD1, CTSW, ARHGDIB, GIMAP4, ISG15, STAT1, KLF6
PC_ 2
Positive: RPL39, TPT1, RPS10, MIF, S100A10, LGALS1, S100A11, S100A6, TMSB10, PGK1
CD7, PLP2, TPI1, GPR183, PCSK1N, TMSB4X, LGALS3, ALDOA, LTB, LSR
RPL21P16, RPS10-NUDT3, CCR7, SPOCK2, CAPG, PLEKHA7, RPL36A, ENSG00000255508, UNC119, RPL21
Negative: CD38, ZEB2, ATXN1, CDK6, ELMO1, HLA-DRB5, CBLB, RABGAP1L, HLA-DRA, FUT8
LYST, PARP14, MIR4435-2HG, NEAT1, PRKCB, SIK3, PAM, RNF19A, RUNX2, ANKRD44
MBNL1, RBPJ, TIAM2, ENSG00000288632, ENTPD1, FYN, UTRN, SFMBT2, RNF213, TOX
PC_ 3
Positive: MALAT1, BTG1, PNRC1, BNIP3L, REL, IL2, PPIAP22, TRIM22, PKIA, RPL21P75
RPL21P119, NINJ1, SLC2A3, LINC01619, RPL23AP42, PCED1B-AS1, ZBTB20, RPL21P16, MAML2, YPEL3
RPL7P9, RPL15P3, HNRNPA1P7, ENSG00000218426, BICDL1, RPL9P7, KLF2, IL7R, ENSG00000280441, P4HA1
Negative: ACTG1, TUBA1B, ACTB, TUBB, PSME2, CORO1A, MCM7, NME1-NME2, MCM5, GINS2
CALR, SLC25A5, DCTPP1, CCT6A, SRM, PA2G4, PRMT1, RAN, HSPE1, HSPD1
RANBP1, DYNLL1, GSTP1, SH3BGRL3, PCNA, DUT, LDHB, NHP2, MCM2, NPM1
PC_ 4
Positive: SLC4A10, IFNG-AS1, PHLDA1, GNA15, IFNG, IL4I1, SLAMF1, HIF1A, SATB1, IL2RA
LGALS3, SLC2A3, RORA, PGK1, CD40LG, CYTIP, NFAT5, FURIN, BHLHE40, IL12RB2
DPP4, SLC16A3, ANKRD28, SYTL3, TNFAIP3, PFKFB3, NFKBIA, KLRB1, GPR183, IL2RB
Negative: IFITM1, RPS4Y1, ISG15, CTSW, CD8B, HLA-DPB1, LY6E, HLA-DPA1, HLA-DRB5, HLA-DRA
HLA-DRB1, PLAAT4, OASL, COTL1, CD27, KLRC4-KLRK1, SP110, IFIT3, IFIT1, C12orf57
HLA-DQA1, CD38, CCDC85B, OAS1, LIMD2, CD52, RPS10-NUDT3, IFIT2, HLA-DMA, ISG20
PC_ 5
Positive: NOLC1, DDX21, NCL, HSPD1, HNRNPAB, EIF4A1, NOP16, FKBP4, EBNA1BP2, SRM
NPM1, TOMM40, C1QBP, SLC4A10, RPS17, PRMT1, NOP56, ODC1, ENSG00000249492, KCNQ5
CCT6A, SATB1, TRAP1, RPL17, DCTPP1, PA2G4, IFNG-AS1, CCT2, HSP90AB1, FEZ1
Negative: SH3BGRL3, C4orf3, S100A10, TPI1, ALDOA, BNIP3L, SIT1, LBH, GPI, IGFLR1
SRGN, CRIP1, GNA15, BNIP3, JUN, LGALS1, MIF, LSP1, CTSB, LDHA
LAG3, APOBEC3G, CD82, PGK1, NMB, TMSB4X, SLC16A3, BTG1, PGAM1, CD4
ElbowPlot(merged.18279, ndims = 200, reduction = "pca")2.5 Plot unintegrated UMAP
merged.18279 <- RunUMAP(merged.18279,
reduction = "pca",
dims = 1:50,
reduction.name = "umap.unintegrated")Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
15:24:48 UMAP embedding parameters a = 0.9922 b = 1.112
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by 'BiocGenerics'
15:24:48 Read 12057 rows and found 50 numeric columns
15:24:48 Using Annoy for neighbor search, n_neighbors = 30
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by 'BiocGenerics'
15:24:48 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
15:24:52 Writing NN index file to temp file /tmp/RtmpzlkBZG/file37c5d823da5d58
15:24:52 Searching Annoy index using 1 thread, search_k = 3000
15:24:56 Annoy recall = 100%
15:24:57 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
15:25:00 Initializing from normalized Laplacian + noise (using RSpectra)
15:25:00 Commencing optimization for 200 epochs, with 532572 positive edges
15:25:09 Optimization finished
DimPlot(merged.18279, reduction = "umap.unintegrated", group.by = c("Site", "Patient", "Timepoint"))2.6 Call preliminary clusters for the purposes of doublet discrimination
merged.18279 <- FindNeighbors(merged.18279, dims = 1:50, reduction = "pca")Computing nearest neighbor graph
Computing SNN
merged.18279 <- FindClusters(merged.18279,
resolution = 0.25,
algorithm = 2)Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 12057
Number of edges: 457888
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9125
Number of communities: 8
Elapsed time: 1 seconds
DimPlot(merged.18279, reduction = "umap.unintegrated", label = T)2.7 Identify and remove doublets
This uses raw counts as input
2.7.1 Combine RNA layers
merged.18279[['RNA']] <- JoinLayers(merged.18279[['RNA']])2.7.2 Convert seurat to sce and check colData
merged.18279.sce <- as.SingleCellExperiment(merged.18279, assay = "RNA")
merged.18279.sceclass: SingleCellExperiment
dim: 61217 12057
metadata(0):
assays(2): counts logcounts
rownames(61217): 5-8S-rRNA 5S-rRNA ... ZZEF1 ZZZ3
rowData names(0):
colnames(12057): P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT ...
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA
colData names(15): orig.ident nCount_RNA ... seurat_clusters ident
reducedDimNames(2): PCA UMAP.UNINTEGRATED
mainExpName: RNA
altExpNames(0):
colData(merged.18279.sce)DataFrame with 12057 rows and 15 columns
orig.ident nCount_RNA
<character> <numeric>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG SeuratProject 20129.8
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT SeuratProject 19535.9
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC SeuratProject 19411.9
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT SeuratProject 20429.8
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC SeuratProject 19299.9
... ... ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT SeuratProject 2048
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT SeuratProject 1499
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT SeuratProject 2483
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC SeuratProject 1177
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA SeuratProject 1471
nFeature_RNA percent.mt
<integer> <numeric>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG 4947 3.86697
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT 4384 2.45138
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC 4735 2.01548
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT 4958 2.64150
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC 4227 4.15930
... ... ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT 965 4.736328
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT 778 1.667779
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT 1214 0.443012
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC 670 2.378929
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA 900 4.979606
miQC.probability miQC.keep
<numeric> <character>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG 0.0756333 keep
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT 0.0434632 keep
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC 0.0573548 keep
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT 0.0424527 keep
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC 0.0975392 keep
... ... ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT 0.05033876 keep
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT 0.02237223 keep
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT 0.27941627 keep
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC 0.00989129 keep
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA 0.07598866 keep
Patient Site
<character> <character>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG P108 IFNg
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT P108 IFNg
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC P108 IFNg
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT P108 IFNg
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC P108 IFNg
... ... ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT P104 IFNg
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT P104 IFNg
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT P104 IFNg
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC P104 IFNg
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA P104 IFNg
Timepoint IpiCohort
<character> <character>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG PostVax 5mgIpi
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT PostVax 5mgIpi
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC PostVax 5mgIpi
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT PostVax 5mgIpi
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC PostVax 5mgIpi
... ... ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT PostVax 2.5mgIpi
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT PostVax 2.5mgIpi
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT PostVax 2.5mgIpi
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC PostVax 2.5mgIpi
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA PostVax 2.5mgIpi
Assay
<character>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG RNA
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT RNA
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC RNA
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT RNA
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC RNA
... ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT RNA
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT RNA
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT RNA
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC RNA
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA RNA
Sample
<character>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG P108_IFNg_PostVax_5m..
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT P108_IFNg_PostVax_5m..
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC P108_IFNg_PostVax_5m..
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT P108_IFNg_PostVax_5m..
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC P108_IFNg_PostVax_5m..
... ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT P104_IFNg_PostVax_2...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT P104_IFNg_PostVax_2...
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT P104_IFNg_PostVax_2...
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC P104_IFNg_PostVax_2...
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA P104_IFNg_PostVax_2...
RNA_snn_res.0.25
<factor>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG 0
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT 0
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC 0
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT 0
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC 0
... ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT 5
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT 1
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT 1
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC 1
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA 1
seurat_clusters ident
<factor> <factor>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG 0 0
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT 0 0
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC 0 0
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT 0 0
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC 0 0
... ... ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT 5 5
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT 1 1
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT 1 1
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC 1 1
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA 1 1
2.7.3 Run scDblFinder
Set the dbr.sd very high to better allow thresholds to be set based on misclassification rates per sample
merged.18279.sce <- scDblFinder(merged.18279.sce,
samples = "Sample",
dbr.sd = 1,
clusters = "seurat_clusters")2.7.4 Inspect results
# Look at the classes
table(merged.18279.sce$seurat_clusters, merged.18279.sce$scDblFinder.class)
singlet doublet
0 2911 196
1 2220 148
2 2095 258
3 1249 167
4 1212 73
5 751 11
6 474 27
7 253 12
table(merged.18279.sce$Sample, merged.18279.sce$scDblFinder.class)
singlet doublet
P101_IFNg_PostVax_2.5mgIpi_RNA 2953 247
P103_IFNg_PostVax_2.5mgIpi_RNA 1099 172
P104_IFNg_PostVax_2.5mgIpi_RNA 3113 190
P108_IFNg_PostVax_5mgIpi_RNA 4000 283
# Look at the scores
summary(merged.18279.sce$scDblFinder.score) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000525 0.0023453 0.0098138 0.1059070 0.0527471 0.9999471
2.7.5 Save doublet classifications into main Seurat object
merged.18279$doublet_classification <- merged.18279.sce$scDblFinder.class2.7.6 Count singlets and doublets
table(merged.18279$doublet_classification)
singlet doublet
11165 892
table(merged.18279$doublet_classification, merged.18279$seurat_clusters)
0 1 2 3 4 5 6 7
singlet 2911 2220 2095 1249 1212 751 474 253
doublet 196 148 258 167 73 11 27 12
2.7.7 Plot singlets/doublets in UMAP space
DimPlot(merged.18279,reduction = "umap.unintegrated", group.by = "doublet_classification")2.7.8 Subset object to remove doublets and count remaining cells
merged.18279.singlets <- subset(merged.18279, doublet_classification %in% c("singlet"))
dim(merged.18279.singlets)[1] 61217 11165
# Count remaining cells per initial cluster
table(merged.18279.singlets$seurat_clusters)
0 1 2 3 4 5 6 7
2911 2220 2095 1249 1212 751 474 253
2.8 Remove cells with very high nCount_RNA values, set other final QC filters
merged.18279.singlets <- subset(merged.18279.singlets,
subset = nCount_RNA < 40000 & nCount_RNA > 3000 & nFeature_RNA > 1500)
dim(merged.18279.singlets)[1] 61217 10320
2.9 Re-compute PCA
Re-scale data now that it has been subset
merged.18279.singlets[["RNA"]] <- split(merged.18279.singlets[["RNA"]], f = merged.18279.singlets$Sample)
merged.18279.singlets <- FindVariableFeatures(merged.18279.singlets,
assay = "RNA",
layer = "data",
selection.method = "vst",
nfeatures = 5000)
merged.18279.singlets <- ScaleData(merged.18279.singlets, vars.to.regress = "percent.mt")
merged.18279.singlets <- RunPCA(merged.18279.singlets, npcs = 200, verbose = TRUE)2.10 Run Harmony integration
merged.18279.singlets <- IntegrateLayers(merged.18279.singlets,
method = HarmonyIntegration,
orig.reduction = "pca",
new.reduction = "integrated.harmony"
)Warning: HarmonyMatrix is deprecated and will be removed in the future from the
API in the future
Warning: Warning: The parameters do_pca and npcs are deprecated. They will be ignored for this function call and please remove parameters do_pca and npcs and pass to harmony cell_embeddings directly.
This warning is displayed once per session.
Warning: Warning: The parameter tau is deprecated. It will be ignored for this function call and please remove parameter tau in future function calls. Advanced users can set value of parameter tau by using parameter .options and function harmony_options().
This warning is displayed once per session.
Warning: Warning: The parameter block.size is deprecated. It will be ignored for this function call and please remove parameter block.size in future function calls. Advanced users can set value of parameter block.size by using parameter .options and function harmony_options().
This warning is displayed once per session.
Warning: Warning: The parameter max.iter.harmony is replaced with parameter max_iter. It will be ignored for this function call and please use parameter max_iter in future function calls.
This warning is displayed once per session.
Warning: Warning: The parameter max.iter.cluster is deprecated. It will be ignored for this function call and please remove parameter max.iter.cluster in future function calls. Advanced users can set value of parameter max.iter.cluster by using parameter .options and function harmony_options().
This warning is displayed once per session.
Warning: Warning: The parameter epsilon.cluster is deprecated. It will be ignored for this function call and please remove parameter epsilon.cluster in future function calls. Advanced users can set value of parameter epsilon.cluster by using parameter .options and function harmony_options().
This warning is displayed once per session.
Warning: Warning: The parameter epsilon.harmony is deprecated. It will be ignored for this function call and please remove parameter epsilon.harmony in future function calls. If users want to control if harmony would stop early or not, use parameter early_stop. Advanced users can set value of parameter epsilon.harmony by using parameter .options and function harmony_options().
This warning is displayed once per session.
Transposing data matrix
Using automatic lambda estimation
Initializing state using k-means centroids initialization
Harmony 1/10
Harmony 2/10
Harmony converged after 2 iterations
2.11 Identify final clusters after integration
merged.18279.singlets <- FindNeighbors(merged.18279.singlets, dims = 1:30, reduction = "integrated.harmony")Computing nearest neighbor graph
Computing SNN
merged.18279.singlets <- FindClusters(merged.18279.singlets,
resolution = 0.7,
algorithm = 2)Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 10320
Number of edges: 352714
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8127
Number of communities: 11
Elapsed time: 1 seconds
merged.18279.singlets <- RunUMAP(merged.18279.singlets,
reduction = "integrated.harmony",
dims = 1:30,
reduction.name = "umap.harmony")15:36:03 UMAP embedding parameters a = 0.9922 b = 1.112
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by 'BiocGenerics'
15:36:03 Read 10320 rows and found 30 numeric columns
15:36:03 Using Annoy for neighbor search, n_neighbors = 30
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by 'BiocGenerics'
15:36:03 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
15:36:06 Writing NN index file to temp file /tmp/RtmpzlkBZG/file37c5d8437d41b
15:36:06 Searching Annoy index using 1 thread, search_k = 3000
15:36:10 Annoy recall = 100%
15:36:11 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
15:36:13 Initializing from normalized Laplacian + noise (using RSpectra)
15:36:13 Commencing optimization for 200 epochs, with 447378 positive edges
15:36:21 Optimization finished
table(merged.18279.singlets$seurat_clusters)
0 1 2 3 4 5 6 7 8 9 10
2023 1984 1845 1165 1093 778 470 373 355 130 104
2.12 Plot clusters
DimPlot(merged.18279.singlets,reduction = "umap.harmony", label = TRUE)DimPlot(merged.18279.singlets,reduction = "umap.harmony", group.by = "Patient")DimPlot(merged.18279.singlets,reduction = "umap.harmony", group.by = "Site")DimPlot(merged.18279.singlets,reduction = "umap.harmony", group.by = "Timepoint")DimPlot(merged.18279.singlets,reduction = "umap.harmony", group.by = "IpiCohort")DimPlot(merged.18279.singlets,reduction = "umap.harmony", group.by = "Sample") + NoLegend()2.13 Save clustered object
saveRDS(merged.18279.singlets,"IFNg_scRNA_Part2.rds")2.14 Get session info
sessionInfo()R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Rocky Linux 8.10 (Green Obsidian)
Matrix products: default
BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so; LAPACK version 3.9.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: America/New_York
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] scDblFinder_1.14.0 SingleCellExperiment_1.22.0
[3] SummarizedExperiment_1.30.2 Biobase_2.60.0
[5] GenomicRanges_1.52.1 GenomeInfoDb_1.36.4
[7] IRanges_2.34.1 S4Vectors_0.38.2
[9] BiocGenerics_0.46.0 MatrixGenerics_1.12.3
[11] matrixStats_1.2.0 harmony_1.2.0
[13] Rcpp_1.0.12 alevinQC_1.16.1
[15] patchwork_1.3.0 Matrix_1.6-4
[17] fishpond_2.6.2 SeuratWrappers_0.3.19
[19] miQC_1.8.0 lubridate_1.9.3
[21] forcats_1.0.0 stringr_1.5.1
[23] dplyr_1.1.4 purrr_1.0.2
[25] readr_2.1.5 tidyr_1.3.1
[27] tibble_3.2.1 ggplot2_3.4.4
[29] tidyverse_2.0.0 Seurat_5.1.0
[31] SeuratObject_5.0.2 sp_2.1-3
[33] sctransform_0.4.1 data.table_1.15.0
loaded via a namespace (and not attached):
[1] spatstat.sparse_3.0-3 bitops_1.0-7
[3] httr_1.4.7 RColorBrewer_1.1-3
[5] tools_4.3.1 utf8_1.2.4
[7] R6_2.5.1 DT_0.31
[9] lazyeval_0.2.2 uwot_0.1.16
[11] withr_3.0.0 GGally_2.2.1
[13] gridExtra_2.3 progressr_0.14.0
[15] cli_3.6.2 spatstat.explore_3.2-6
[17] fastDummies_1.7.3 labeling_0.4.3
[19] spatstat.data_3.0-4 ggridges_0.5.6
[21] pbapply_1.7-2 Rsamtools_2.16.0
[23] R.utils_2.12.3 scater_1.28.0
[25] parallelly_1.37.0 limma_3.56.2
[27] generics_0.1.3 BiocIO_1.10.0
[29] gtools_3.9.5 ica_1.0-3
[31] spatstat.random_3.2-2 ggbeeswarm_0.7.2
[33] fansi_1.0.6 abind_1.4-5
[35] R.methodsS3_1.8.2 lifecycle_1.0.4
[37] yaml_2.3.8 edgeR_3.42.4
[39] Rtsne_0.17 grid_4.3.1
[41] promises_1.2.1 dqrng_0.3.2
[43] crayon_1.5.2 shinydashboard_0.7.2
[45] miniUI_0.1.1.1 lattice_0.22-5
[47] beachmat_2.16.0 cowplot_1.1.3
[49] pillar_1.9.0 knitr_1.45
[51] metapod_1.8.0 rjson_0.2.21
[53] xgboost_1.7.7.1 future.apply_1.11.1
[55] codetools_0.2-19 leiden_0.4.3.1
[57] glue_1.7.0 remotes_2.4.2.1
[59] vctrs_0.6.5 png_0.1-8
[61] spam_2.10-0 gtable_0.3.4
[63] xfun_0.42 S4Arrays_1.2.0
[65] mime_0.12 survival_3.5-8
[67] statmod_1.5.0 bluster_1.10.0
[69] ellipsis_0.3.2 fitdistrplus_1.1-11
[71] ROCR_1.0-11 nlme_3.1-164
[73] RcppAnnoy_0.0.22 irlba_2.3.5.1
[75] vipor_0.4.7 KernSmooth_2.23-22
[77] colorspace_2.1-0 nnet_7.3-19
[79] tidyselect_1.2.0 compiler_4.3.1
[81] BiocNeighbors_1.18.0 DelayedArray_0.26.7
[83] plotly_4.10.4 rtracklayer_1.60.1
[85] scales_1.3.0 lmtest_0.9-40
[87] digest_0.6.34 goftest_1.2-3
[89] spatstat.utils_3.0-4 rmarkdown_2.25
[91] RhpcBLASctl_0.23-42 XVector_0.40.0
[93] htmltools_0.5.7 pkgconfig_2.0.3
[95] sparseMatrixStats_1.12.2 fastmap_1.1.1
[97] rlang_1.1.3 htmlwidgets_1.6.4
[99] shiny_1.8.0 DelayedMatrixStats_1.22.6
[101] farver_2.1.1 zoo_1.8-12
[103] jsonlite_1.8.8 BiocParallel_1.34.2
[105] R.oo_1.26.0 BiocSingular_1.16.0
[107] RCurl_1.98-1.14 magrittr_2.0.3
[109] modeltools_0.2-23 scuttle_1.10.3
[111] GenomeInfoDbData_1.2.10 dotCall64_1.1-1
[113] munsell_0.5.0 viridis_0.6.5
[115] reticulate_1.35.0 stringi_1.8.3
[117] zlibbioc_1.46.0 MASS_7.3-60.0.1
[119] plyr_1.8.9 flexmix_2.3-19
[121] ggstats_0.5.1 parallel_4.3.1
[123] listenv_0.9.1 ggrepel_0.9.5
[125] deldir_2.0-2 Biostrings_2.68.1
[127] splines_4.3.1 tensor_1.5
[129] hms_1.1.3 locfit_1.5-9.8
[131] igraph_2.0.2 spatstat.geom_3.2-8
[133] RcppHNSW_0.6.0 reshape2_1.4.4
[135] ScaledMatrix_1.8.1 XML_3.99-0.16.1
[137] evaluate_0.23 scran_1.28.2
[139] BiocManager_1.30.22 tzdb_0.4.0
[141] httpuv_1.6.14 RANN_2.6.1
[143] polyclip_1.10-6 future_1.33.1
[145] scattermore_1.2 rsvd_1.0.5
[147] xtable_1.8-4 restfulr_0.0.15
[149] svMisc_1.2.3 RSpectra_0.16-1
[151] later_1.3.2 viridisLite_0.4.2
[153] beeswarm_0.4.0 GenomicAlignments_1.36.0
[155] tximport_1.28.0 cluster_2.1.6
[157] timechange_0.3.0 globals_0.16.2